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Abstract 

We describe algorithms for finding harmonie cochains, an essential in- 
grédient for solving elliptic partial differential équations using finite élé- 
ment or discrète exterior calculus. Harmonie cochains are also useful in 
computational topology and computer graphies. We focus on finding har- 
monie cochains cohomologous to a given cocycle. Amongst other things, 
this allows for localization near topological features of interest. We dérive a 
weighted least squares method by proving a discrète Hodge-deRham the- 
orem on the isomorphism between the space of harmonie cochains and 
cohomology. The solution obtained either satisfies the Whitney form fi- 
nite élément exterior calculus équations or the discrète exterior calculus 
équations for harmonie cochains, depending on the discrète Hodge star 
used. 
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1 Introduction 

We discuss methods for finding simplicial harmonie cochains - approximations 
of harmonie forms on simplicial meshes. In particular, we want to find the har- 
monie cochain cohomologous to a given cocycle. That is, given a cocycle ù), we 
want a harmonie cochain h such that h = ù) + dafoY some a. We either solve an 
eigenvector problem foUowed by post processing or use a weighted least squares 
method. 

Harmonie cochains are used in finite élément solution of elliptic partial dif- 
ferential équations like the Poisson's équation ApU = f. See for instance [2]. 
They are also useful in computer graphies for design of vector fields, since they 

* This is a much shorter incarnation of version 6 of this paper which is available on arXiv as [U] . 
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can provide a background on which vortices, sources and sinks may be super- 
imposed [9] . In computer graphies they are also useful for finding conformai 
parameterization for texture mapping and other applications [10] . 

We prove an easy discrète version of the Hodge-deRham isomorphism the- 
orem. This leads to a weighted least squares based method which is the main 
contribution of this paper. The linear System is an obvions one and can be de- 
rived also from the gradient part of Hodge décomposition or in other ways. The 
two other methods we describe are based on finding eigenvectors foUowed by 
post processing. The least squares method solves the mixed finite élément ex- 
terior calculus équations for harmonie cochains given in [2, Lemma 3.10]. (This 
is a resuit of Demlow and Hirani, and the proof can be found in [11].) For each 
of the harmonie cochain methods considered, the choice of the Hodge star op- 
erator (Whitney or primal-dual) can be made, leading to two variations of each 
method. 

Other methods are those by Gu and Yau [10] , and Desbrun et al. [7] . Both of 
thèse have some numerical disadvantages especially when Whitney Hodge star 
is used instead of the diagonal primal-dual Hodge star of discrète exterior cal- 
culus. (The Whitney Hodge star is needed for gênerai simplicial meshes, and for 
the lowest order finite élément exterior calculus.) In cases such as 2-dimensional 
cochains in tetrahedral meshes, the Desbrun et al. method does more work than 
is necessary for forming the linear System, no matter which Hodge star is used. 

2 Prelîminaries 

]V[ost of the needed background information on algebraic topology and exterior 
calculus can be found in an earlier longer version of this paper which is still 
available on arXiv [11]. We use two types of discretizations of exterior calculus - 
discrète exterior calculus, and finite élément exterior calculus. In finite élément 
exterior calculus, we only consider the version that uses Whitney forms. 

We first recall the smooth Hodge-deRham theorem on the isomorphism be- 
tween cohomology and harmonie forms (ker A) or harmonie fields (ker d n ker5). 
(This material is based on [14]). The space of harmonie p-dimensional fields on 
a manifold M is denoted ^^(M). For a elosed manifold (i.e., compact mani- 
fold without boundary), harmonie forms and harmonie fields are the same, i.e., 
ker A = kerdnker5. However, in the case of compact manifolds with boundary 
ÔM, which we will refer to as (9-manifolds, one only has that ker d n ker 5 e ker A 
and there can exist harmonie forms which are not harmonie fields [6] . 

One of the striking properties of harmonie forms or fields is the link they 
yield between topology and analysis or geometry. For elosed manifolds there 
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is an isomorphism between real cohomology and the space of harmonie forms. 
For compact 5-manifolds however, even the space of harmonie fields is infinité 
dimensional due to the possibility of specifj^ing boundary conditions. An iso- 
morphism with cohomology can be obtained by restricting harmonie fields by 
specifying certain boundary conditions. 

The tangential component of a p-form o) is denoted tœ and its value is the 
value of (X) on the tangential (to dM) components of its vector field arguments. 
Thenthe norma/ component ofo^isno) = ù)\qm-^(^' See [14, page 27] or [1, page 
540] . Thèse can also be defined using the puUback via the inclusion map of the 
boundary into the manifold. A differential form ù) is said to satisfy the Neumann 
or absolute boundary conditions if it has zéro normal component (no) = 0), and 
the Dirichlet or relative boundary conditions if it has zéro tangential component 
{X(x) = 0). Let ^^(M) and (M) be harmonie fields satisfying the Neumann or 
Dirichlet boundary conditions, respectively. Then one has: 

Theorem (Hodge-deRham Isomorphism [14]). If M is a closed manifold, then 
Hf^{M'M) =J^f^{M) =kerAp, and if it is a compact d -manifold then HP{M;U) = 
^^(M) andHP{M,dM;m=^D^M). 

The space H^iM'M) is the (absolute) real p-cohomology vector space of M, 
and H^{M,dM;U) is the relative real p-cohomology vector space of M, relative 
to its boundary. For ô-manifolds, we will only consider harmonie fields satisfy- 
ing Neumann conditions. This is because the least squares method is based on 
a weak form of the Laplace-deRham operator, and in that framework the Neu- 
mann conditions are automatic, that is they do not have to be enforced explicitly. 
For manifold complexes with boundary we will use harmonie cochains synony- 
mously with harmonie Neumann cochains. 

3 Eîgenvector Methods 

Cohomologous harmonie cochains can be computed by first Computing a har- 
monie cochain basis foUowed by some post processing. Such a basis can be 
obtained as eigenvectors of the zéro eigenvalue of a discrète A^. The problem 
of finding eigenvectors can be formulated (in the terminology of finite élément 
methods) using a weak mixed or weak direct method. While nothing is pub- 
lished about the eigenvector method, the weak mixed method was the one used 
by Arnold et al. [2] in one of their examples. 

Let A = d5 + (5d be the smooth Laplace-deRham operator on some mani- 
fold M. Then the direct eigenvalue problem is to find a nonzero differential 
form u and a real scalar A such that Au = Xu. The formai dérivation of the 
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weak direct method goes like this: start by posing the problem of finding a u 
such that (A u^v) = X{u, v) for ail v, the inner products being those on forms. 
Then using the formula for the Laplace-deRham operator, and assuming ap- 
propriate boundary conditions (which implies adjointness of d and ô) this is 
équivalent to finding a u such that (d u,dv) + {Ô Uyôv) = A( v) for ail v. If M 
is replaced by its simplicial complex approximation (which we will also refer to 
as M) then the discretization yields the linear System /\pU = X^pUy where now 

:= dp ^p+idp + C-l)^^"^^^""^"^^^ *pdp_i *~iid^_j *p is the discrète Laplace- 
deRham operator [11] and w is a p-cochain. Here *p is the mass matrixfor Whit- 
ney p-forms or the primal-dual discrète Hodge star. The harmonie cochains 
are thus the solutions corresponding to the zéro eigenvalue for this generalized 
eigenvalue problem. 

For the weak mixed eigenvector method, consider the linear System for the 
unknowns a and w. 

(o-,T)-(dp_iT, l/) = 0, 
(dp_i a, v) + (dp w, dp [;) = , 

for ail T and v. Then (a, u) is a solution if and only if a = and w is a harmonie 
p-form [2, Lemma 3.10]. We discretize thèse équations and obtain the System 
matrix 

*pûp_i d^*p+idp 

whose eigenvectors corresponding to the zéro eigenvalue we seek. 

Figure 1 shows results of the eigenvector calculations.The eigenvector meth- 
ods will often suffice, if ail that is needed is some harmonie basis, which may be 
the common case in finite élément exterior calculus. Applications like vector 
field design in computer graphies may require more control over the process, 
namely the satisfaction of the cohomology constraint. 

3. 1 Projection based methods 

If a harmonie cochain basis is available, then orthogonal projection to the har- 
monies can be used to obtain a harmonie cochain h cohomologous to a given 
cocycle œ, (This method was suggested to us by Ari Stern.) In contrast, the least 
squares method discussed in Section 4 finds a cohomologous harmonie cochain 
without requiring any precomputation of a harmonie basis. IMoreover, the pro- 
jection method does not find the potential of the gradient part. If that is needed, 
then the least squares method équation (4) has to be solved anyway. 
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Figure 1: Harmonie cochains produced by the mixed eigenvector method. The torus 
has a two-dimensional space of harmonie cochains and the four-holed dise has a four- 
dimensional space of harmonie Neumann fields. 
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Let Hhea matrix whose columns form a harmonie p-cochain basis. Given 
a nontrivial cocycle o), we seek the harmonie cochain h such that h = co + dafor 
some a. (Thus we are interested in a Hodge deeomposition of co. Note that the 
Hodge décomposition of an arbitrary co would be d a + 5 ^ + /z, but since the given 
ù) is 3. nontrivial cocycle, it has no curl part.) Since imd is orthogonal to every 
column hi ofH, wehave that [ùj + da.hi) = [co, ht) - [Y.j ajhj.hi) for ail where 
X; cijhj is the h that we seek. (The inner product above is the p-cochain inner 
product.) Writing a for the vector of unknown coefficients aj, we can express 
the last equality as the linear System ^ Ha = * o). After solving this for the 
unknowns a, the vector Hais the desired h. This is the normal équation for a 
weighted least squares problem (a différent System from theh one in Section 4). 
The matrix of the linear System is of order of the p-Betti number and the cost of 
this projection will be dominated by the matrix vector multiplications needed in 
forming ^ Hiî the Betti number is small. If the columns of H are orthonormal 
in the * inner product then no linear solve is required. 

3.2 Paîrîng with homology basis 

For vector field design in computer graphies or in physical applications, the 
usual cases are dimension 2 with 1-cochains and dimension 3 with 1-cochains 
or 2-cochains. In the latter case, only solid handles and cavities are relevant 
since gênerai 3-manifolds are typically not used in such applications. In ail thèse 
cases, it makes sensé to talk of homology basis éléments corresponding to topo- 
logical features. Thèse can be used by pairing with cohomology to find coho- 
mologous harmonie cochains. (This method was suggested to us by Douglas 
Arnold.) For this one needs an explicit isomorphism H^[M) = //p(M)*, where 
Hp (M) * is the vector space dual of real-valued homology. Note that this method 
requires not only the entire harmonie cochain basis but also a fuU homology 
basis. 

Given [œ] eH^iM). define the map cp \ {M) ^ Hp{M)'' by ^[do] [z] \=ù){z) 
for any [z] g Hp (M) . This map is well-defined: given other représentatives œ+Aa 
and z + dyy one has [ù) + Aa)[z + dy) = ù){z) + da{z) + ù){dy) + da{dy) = ù){z). To 
prove that cp is an isomorphism, it is enough to show that it is injective. That 
is, we would like to show that for any [oj] e H^[M), ù)[z)=0 for ail [z] e Hp[M) 
implies that [ù)] = 0. This is équivalent to showing that if is a représentative 
of an élément of H^{M), ù){z) = for ail nontrivial cycles z implies that o) is 
exact. Since ù) is nontrivial, ù) = da + hfor some a and harmonie cochain h. To 
show ù) is exact is the same as showing h = 0. Thus we have to show that given 
a harmonie cochain h, hiz) = for ail nontrivial cycles z implies that h = 0. We 
now show this for the case of n = 2 (dim of M) and p = l. 
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Theorem 1. LetM be a surface simplicial complex. Then cp : (M) ^ Hi (M) * is 
injective (hence an isomorphism). 

Proof. It is enough to consider a homology basis of nontrivial cycles. Suppose 
M has b holes and g handles. Consider a homology basis corresponding to the 
holes, handles and tunnels. That is, let zi , . . . , be cycles corresponding to 
- 1 of the b holes (the remaning one hole is considered the outer boundary), 
jUi , . . . , jUg be handle cycles corresponding to the g handles, Ai , . . . , Ag be tunnel 
cycles corresponding to the g handles. (Handle cycles are like longitudes on a 
torus and tunnel cycles are like latitudes on a torus.) Let (Oz be the collection 
of - 1 nontrivial cocycles corresponding to the hole cycles in, and let o)^ and 
ù)x be similarly defined. Each such cocycle cois s. "picket fence" (see Figure 2). 
Either two edges of a triangle carry a part of co or none. The hole cocycles join 
the boundary of a hole to the outer boundary. The handle and tunnel coycles go 
around the handle or tunnel. Such cocycles are obtained by dualizing cycles on 
the dual mesh. By Theorem 2, there is a basis of cohomologous cochains hz for 
ail z, for ail jU, and hx for ail A (cohomologous to the corresponding œ's). 

Now consider a harmonie 1-cochain that évaluâtes to on ail the basis cycles 
above. In terms of the harmonie basis above, 

h = Y.^zhz + Y.^^ih^ + Y.^xhx' (2) 

z 11 X 

Note that hz évaluâtes to nonzero on z and on every other cycle, h^. éval- 
uâtes to nonzero on A/ and on every other cycle, hx^ évaluâtes to nonzero on 
jii and on every other cycle. hz[z) = cOziz) = cOziBz), where Bz is the bound- 
ary of the corresponding hole since hz is cohomologous to cOz and z is homolo- 
gous to Bz. But cOziBz) = ±1 (or whatever value was picked for edges). Likewise, 
hziz') = o)z{z') = (J^ziBz') = 0, for z' ^ z since ù)z takes value on edges of Bz'. 
Similarly for hz on other types of cycles, and for the other harmonie basis élé- 
ments. Thus, the coefficients in (2) are ail zéro. □ 

If M is the closure of a connected open subset of and the topological fea- 
tures of interest are cavities and solid handles then a resuit similar to the above 
one can be shown. Now let // be a matrix whose columns form a basis of har- 
monie p-cochains and B a matrix whose columns form a homology basis corre- 
sponding to topological features in the sensé described in the proof above. Then 
{B^ H)~^ contains the harmonie cochains cohomologous to the topological fea- 
tures. 
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4 Least Squares Method 



In what foUows, M will be a simplicial manifold complex, with or without bound- 
ary. Ail références to A are to the discrète Laplace-deRham operators [11]. For 
a closed manifold, one way to show the Hodge-deRham isomorphism theorem 
of Section 2 for the smooth case is to use a variational approach [12, Theorem 
2.2.1]. One shows that in each cohomology class there is exactly one harmonie 
form and it is the one with the smallest norm. The norm used is the norm 
induced from the inner product of differential forms. Inspired by this, we for- 
mulate a simple discrète version of this theorem. This is donc for harmonie 
cochains in the case of manifold simplicial complexes without boundary, and 
for harmonie Neumann cochains in the case with boundary. First we dérive the 
necessary stationarity conditions in the discrète case. For coeC^ s.t. dpù) = 0, we 
consider the optimization problem min^j^^cp i + dp_i a, co + dp-i a)cp, where 
the (•, •)cp is the inner product on p-cochains [3] . Writing this in matrix notation, 
we want to find the minimizer a in the optimization problem 

min (ù) + dij-i a]^ ^ D (ù) + di)-i a] . (3) 

From the stationary condition for the minimizer and using properties of the 
Hodge star matrix, we obtain: 

d^_i *pdp-ia = -d^_i ^pùj. (4) 

This is the normal équation for the weighted least squares problem da--ù). 
Although the above équation is a necessary condition for soMng the optimiza- 
tion problem (3), the matrix d^_^ ^pdp-i may have a nontrivial kernel. In fact in 
the interesting cases it generally will. (For example, for p = 1, the kerdo will 
have dimension equal to the number of connected components in the com- 
plex.) Thus, for a to be a minimizer we need that the Hessian d^_^ ^pdp-i be 
at least positive semidefinite, which is true because of the positive definiteness 
of In this case, a may not be unique, but as we will show next, dp-i a will 
be unique. Note that équation (4) is équivalent to ôpdp-i a = -ôpCO which is 
ôp{ù) + dp-i a) = 0. This should make the connection to ù) + dp-i a being har- 
monie more transparent since we also have that dio) + da) = 0. 

From the above, if [co] e H^iK'M), then it is easy to see that (i) there exists 
a cochain a e C^~^iK; U), not necessarily unique, such that ôp [œ + dp-i a) = 0; 
(n) there is a unique cochain dp-i a satisfying ôp {oj+dp-i a) = ; and (Ui) ôp {ù)+ 
dp-i a) = implies {o) + dp-i a) = 0. To see (i) consider the least squares 
problem dp-i a-co. Let -a be a solution. Some such a always exists because 
least squares problems always have a solution. Note that the norm used in 
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formulating this problem as a residual minimization is the one induced from 
the Hodge star inner product on cochains. Specifically, the inner product ma- 
trix is *p and the least squares problem minimizes [o) + dp-i a)^ ^p{ù) + dp-i a) 
since co - dp-i (-a) = {(0 + dp-i a) is the residual. But from properties of least 
squares [5] the residual {o) + dp-i a) is *p -orthogonal to imdp_i. Thus we have 
that iù) + dp-i a) e imdp-i"'"*^ = kerÔp since ôp is the adjoint of dp-i up to sign 
in the Hodge star inner product on cochains. In (ii) uniqueness of dp-i a fol- 
lows from properties of least squares, and (iii) is obvions since co + dp-i a is also 
closed. Note that unlike in the smooth case, Ôp+i and dp are adjoints of each 
other up to sign only. Specifically, [dp a, 13)^^^^ = (-1)^"^ [a,Ôp+il3)^^ for any 
p-cochain a and {p + l)-cochain jS. From the preceding discussion, we have the 
foUowing elementary but useful theorem: 

Theorem 2 (Discrète Hodge-deRham Isomorphism). There is a unique harmonie 
eoehain in eaeh eohomology elass and it is the one with the smallest norm. Given 
a cocycle co its cohomologous harmonie eoehain isco + da where a is a solution of 
d^*da = -d^ ^ù). 

An alternative dérivation of (4) is to project co to image of d by requiring that 
(da,dT) = {ù),dT) for ail t. Yet another dérivation is the foUowing. Given an 
ùj, to find its Hodge décomposition, one starts with co = da + ô 13 + h, where we 
are seeking a harmonie field or eoehain h and an a and j6. Applying ô to both 
sides yields *"M^*da=*"M^*a), which is the same as (4) up to sign after the 
is cancelled from both sides. Note that the linear System for the j6 part is 
dô 13 = dùj. This has a which cannot be removed by cancellation since this 
is d j6 = dù). In his thesis [4], Bell was motivated by the need to address 

the inverse Hodge star matrix in order to apply algebraic multigrid to the Hodge 
décomposition problem. He proposed replacing the Hodge stars by identity and 
solving the above Systems starting with random cochains until one has obtained 
a eohomology basis. (He did not prove that the procédure is guaranteed to pro- 
duce such a basis.) He then showed that choosing the basis éléments as ù) and 
solving (4) for each one yields a basis of harmonie cochains. In contrast we have 
shown above that each such o) is cohomologous to the corresponding harmonie 
eoehain h individually 

Ail computations in this paper were donc using the Python language with 
SciPy, NumPy, and PyDEC [3] packages. The top two rows of Figure 2 show 
the harmonie cochains cohomologous to given nontrivial cocycles on a torus 
surface. The bottom two rows of Figure 2 show several examples on a planar 
mesh with holes. To single out a particular hole, so that the harmonie eoehain 
proxy vector field will circulate around that hole, one picks a cocycle Connecting 
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Figure 2: Some example computations using the least squares method. Top two rows : 
Cocycles representing a cohomology basis for the torus are shown as thick edges in the 
left figures. Thèse are the œ cocycles of the text. The cocycles have value ±1 on thèse 
edges and on the other edges. The right figures show the harmonie cochains in the cor- 
responding cohomology classes. Bottom two rows : The nontrivial cocycles are marked 
in red. Note that the proxy vector fields circulate around only those holes associated 
with the co cycle, and past others. 
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Figure 3: Three différent cocycles {œ of the text) representing the same cohomology 
class lead to the same harmonie cochain when least squares method is used. 

that boundary to the outer boundary. Connecting two holes results in a har- 
monie cochain that circulâtes about those two holes. For the cochains shown 
in the third row of Figure 2, from left to right, the values of ||A/z|Uj relative to 
||/z|U^ are approximately 4.12 x 10~^\ 7.32 x 10"^^ and 3.02 x 10~^\ respectively. 
Similarly, for the cochains in the bottom row, from left to right, thèse values are 
4.98 X 10~^\ 5.73 x 10"^^ and 5.64 x 10~^\ respectively. Figure 3 shows that the 
least squares method (as expected) finds the same harmonie cochain when very 
différent initial cocycles from the same cohomology class are given as input. If 
the cohomologous cochains are denoted h, h' , h" from left to right, respectively, 
then the différences between them are \\h - h'\\^^ = 1.1 x 10"^"*, \\h - h"\\^^ = 
2.8 X 10-14 and || h' -h"\\^^= 2.2 x 10"!^ 

4. 1 Lînear solvers for the least squares method 

As noted earlier, the matrix d^_^ * p dp_i in (4) is positive semidefinite since dp_i 
will typically have a nontrivial kernel. For example, for p = 1 for a connected do- 
main, the space of constant functions on the domain is in the kernel of do- In 
this case, it is easy to make the System nonsingular (mod out the nontrivial ker- 
nel) by fixing the value at a vertex and adjusting the linear System accordingly. 
For the case of 2-cochains in tetrahedral meshes however, the kernel of di can 
be large. Let M be a three-dimensional manifold simplicial complex. Simple 
linear algebra and elementary topology reveals that the dim(kerdi) > Nq- xi^) 
where Nq is the number of vertices and xiK) is the Euler number (the alter- 
nating sum of Betti numbers at ail dimensions) [13]. For example, for a con- 
nected domain with boundary, we will have dim(kerdi) > number of vertices - 
1 + number of solid handles - number of cavities. By refining the mesh this ker- 
nel dimension can be made arbitrarily large. If a direct solver is to be used for 
solving (4) then one must mod out this potentially large nontrivial kernel. An 
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alternative is to use itérative Krylov solvers as they work well even in the prés- 
ence of a nontrivial kernel and this is the approach we chose in our experiments. 
Specifically, we used a conjugate gradient solver without any preconditioning or 
modifications. Algebraic multigrid is another very efficient alternative whose 
effectiveness for this problem has been shown in [4]. 

4.2 Finding the initial nontrivial cochains 

In this paper we assume that a nontrivial cocycle is given. Our aim here is not 
to give algorithms for finding a cocycle. However, a few words about this are 
in order. An initial nontrivial cocycle in a cohomology class can be found in a 
number of ways. For surfaces, efficient algorithms to do this exist. By a folklore 
theorem, in time linear in the number of simplices, one can find a homology 
basis for the topological dual (e.g., barycentric dual) graph of the triangulation. 
One can then use Poincaré-Lefschetz duality [13] to get a cohomology basis on 
the primai mesh. For a boundaryless manifold simplicial complex, one would 
start with nontrivial cycles on the dual graph. But in case of a manifold with 
boundary, due to Lefschetz duality, one has to start with a nontrivial relative cy- 
cle on the dual mesh, relative to the boundary. One can also start with a random 
cochain and compute the desired nontrivial cocycle using a Hodge décomposi- 
tion with standard inner product [4]. Yet another method is to use the persis- 
tence algorithm [8]. This is usually implemented using coefficients in finite field 
F2 and has cubic (in the number of simplices) complexity. 

5 Comparisons with Other Methods 

The first relevant method to compare with is from the book of Gu and Yau [10] 
and also appears in their earlier work. The formulation is very simple and straight 
forward, but it leads to inefficient methods on gênerai simplicial meshes. This 
method was further simplified by Desbrun et al. [7] who solve a Poisson's-like 
équation at a différent dimension. The resulting linear Systems in both methods 
suffer from numerical and scalability issues for gênerai simplicial meshes. 

Gu and Yau start with a nontrivial cocycle co representing a cohomology class 
in H^{K) and seek a cochain ù) + àa' such that A(a) + daO = 0. This leads to the 
linear System dp_i d^_^ *pdp_i a' = -dp-i d^_^ ^pCO. The présence of 
the inverse Hodge stars in this Systems lead to numerical disadvantages. 

Desbrun et al. [7] solve a différent Poisson's équation {dp-2 ôp-i + ôp dp-i)a^ 
= -Ôpù). A solution a' to the above équation yields anco + da' that is harmonie 
in the sensé of this paper. Of course, if harmonie 1-cochains are being sought, 
then a' is a 0-cochain and Ôq is the operator. Thus the dô term is not présent. 
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However, the dÔ term is superfluous at every dimension as we have shown. Thus 
their linear System has an extra, unnecessary term. This extra term causes nu- 
merical and scalability problems when Whitney Hodge star is used. 

In Figure 4, we compare the sparsity of the least squares and Desbrun et 
al. matrices for finding harmonie 2-cochains on a tetrahedral mesh of a solid 
annulus (a solid bail with an internai cavity). The matrices are shown in Fig- 
ure 4 for both the Whitney and DEC Hodge stars. Both matrices are of the same 
size but the Desbrun et al. matrix is denser. This is very obvions for the Whitney 
Hodge star case (14.3 million vs. 56 thousand nonzeros). However, it is also évi- 
dent in the DEC Hodge star case (94 thousand vs. 40 thousand nonzeros). Here 
the increased density is due to the extra term in the Desbrun et al. System. 

The superior sparsity of the linear System matrix in the least squares method 
leads to improved solution time. To illustrate this, we compare the time taken 
for again finding harmonie 2-cochains on a tetrahedral mesh of a solid annulus. 
For the least squares method, using conjugate gradient method (without pre- 
conditioning), the times are 0.1355 and 0.1181 seconds for the DEC and Whit- 
ney Hodge stars, respectively. For the Desbrun et al. method, thèse times are 
3.510 and 1746 seconds, respectively. We also used a sparse solver in SuperLU for 
Desbrun et al. System and in this case, the times are 0.3171 and 13.05 seconds, 
respectively. (Ail times are averaged over many trials. Also, it may be possible 
to improve the times for both the methods by using preconditioners or spécial 
solvers.) Another least square method is that of Fisher et al. [9] . Comparisons 
with it are in an earlier version of this paper available on arXiv [11]. 

6 Conclusions 

We presented two methods for finding harmonie cochains in the cohomology 
class of a given cocycle - an eigenvector method (using direct or mixed for- 
mulation) foUowed by post processing and a least squares method. The most 
salient feature of the least squares method is in finding a cohomologous har- 
monie cochain without requiring an entire harmonie or homology basis. The 
least squares method is numerically superior and independent of the choice of 
Hodge stars in comparison with the Poisson's équation methods of Ou and Yao, 
and Desbrun et al. In future we plan to develop harmonie cochain methods for 
higher order finite élément exterior calculus analogous to the one for Whitney 
forms. A précise quantification of the efficiency of the least squares method in 
comparison with the eigenvector method for finding a cohomologous harmonie 
basis is another direction to pursue. 
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Least squares matrix Desbmn et al. matrix 




Figure 4: Magnitudes of nonzeros in operators using the Whitney (top row) and DEC 
(bottom row) Hodge stars. The least squares matrix (left column) is sparser than the 
Desbrun et al. matrix (right column) in the case of DEC Hodge star and significantly 
sparser in the case of Whitney Hodge star. This is due to the extra term in the Desbrun 
et al. matrix. The colorbar shows the magnitude of the nonzero components. The two 
matrices are of equal size, and are for finding harmonie 2-cochains on the tetrahedral 
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